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1.0  INTRODUCTION 


In  recent  years  some  important  developments  have  occurred  in  the  theory  of  two- 
dimensional  laminar  boundary- layer  flows.  The  crucial  discovery,  due  to  van 
Dommelen  and  Shen1,  is  that  the  solution  of  the  boundary- layer  equations  with 
boundary  conditions  corresponding  to  a  circular  cylinder  started  impulsively 
from  rest  develops  a  singularity.  This  with  u^  denoting  the  velocity  of  the 
cylinder  and  a  its  radius,  occurs  at  a  time  u^t/a  =  3.0  and  at  an  angular 
distance  9  »  111°  from  the  forward  stagnation  point.  At  this  time  the  position 
of  zero  skin  friction  is  at  9  =  106°  close  to,  but  not  coincident  with,  the 
singularity.  This  discovery  was  made  by  solving  the  governing  boundary- layer 
equations  in  Lagrangian  form;  subsequently,  the  existence  of  the  singularity 
was  confirmed  by  Cebeci  using  the  Eulerian  form  and  by  Cowley  using  the 
method  of  series  truncation.  The  importance  of  this  result  is  two-fold. 

First,  it  is  of  importance  to  the  fundamental  theory  of  high-Reynolds  number 
flows  in  indicating  that  an  interactive  theory  is  necessary  to  understand  the 
evolution  of  the  flow  field  at  a  finite  time  after  the  motion  starts,  rather 
than  only  after  a  long  time.  Secondly,  it  is  clear  that  significant  changes 
will  then  occur  in  the  flow  field  and  these  may  be  related  to  the  phenomenon 
of  vortex  shedding.  The  physical  processes  which  take  place  during  the  oscil¬ 
lation  of  the  angle  of  attack  of  the  airfoil  are  complicated  and  depend  on  a 
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large  number  of  parameters  *  .  For  example,  an  important  characteristic  is 
a  large  vortex  that  Is  formed  near  the  surface  at  some  stage  in  the  cycle  and 
causes  stall  to  occur  shortly  afterwards.  The  occurrence  of  the  vortex  Is 
probably  associated  with  a  breakdown  of  the  unsteady  boundary  layer®. 

The  purposes  of  the  present  research  are  to  determine  the  relationship  between 
unsteady  separation  and  singularities  in  the  solution,  and  to  explore  the 
possibilities  of  removing  this  singularity  by  interaction  of  the  viscous  and 
Invlscid  equations.  In  this  report  we  describe  the  work  that  has  been  accomp¬ 
lished  under  the  contract.  So  far  we  have  examined  the  evolution  of  the 
boundary  layer  near  the  nose  of  an  oscillating  airfoil  and  found  that,  when  the 
reduced  frequency  Is  of  the  same  order  as  In  the  experiments  on  dynamic  stall. 


the  unsteady  boundary  layer  ceases  to  behave  In  a  smooth  manner  just  downstream 
of  separation  and  before  one  cycle  has  been  completed;  as  with  the  Impulsively 
started  circular  cylinder,  this  Irregular  behavior  signals  the  onset  of  a 

singularity  In  the  solution  of  the  boundary- layer  equations.  The  equations  and 
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solution  procedures  used  in  the  present  study  are  described  in  Sections  2  and 
3,  respectively,  and  the  results  are  presented  and  discussed  in  Section  4. 


2.0  EQUATIONS,  INITIAL  AND  BOUNDARY  CONDITIONS 


The  boundary- layer  equations  and  their  boundary  conditions  for  unsteady 
incompressible  laminar  flows  on  oscillating  airfoils  can  be  written  as 


3u  3v 
Is  +  3n  = 


0 


au  +  u  9u  +  v  3u  =  3Ue  +  u  3Ue  +  v  *Z[1 
Tf+uTs  +  vTn  TF"  +  ue  +  v  TT 

3n 


u(s,0,t)  =  0,  v(s,0,t)  =  0,  u(s,ne,t)  =  ug(s,t) 


(2.1) 

(2.2) 

(2.3) 


The  solution  of  boundary- layer  equations  requires  that  the  external  velocity 
distribution  be  specified.  Since  the  present  effort  is  directed  toward  solu¬ 
tions  near  the  leading  edge  of  the  airfoil,  a  local  model  for  the  potential 
flow  has  been  chosen  in  the  place  of  a  full-potential-flow  code.  We  consider 
an  ellipse  with  major  axis  2a  and  minor  axis  2ax(x«l)  at  an  angle  of 
attack  a(t).  The  surface  of  the  body  is  defined  by 


x  =  -a  cose,  y  =  ax  sine  -n  <  e  <  n 


and  with  these  definitions  and  to  a  first  approximation,  the  external  velocity 
around  the  ellipse  can  be  deduced  from  inviscid  flow  theory  to  be 

o  5  + 

u°(s,t)  - - (2.4) 

e  - 7 

/i  +  r 

Here  u^(s,t)  denotes  a  dimensionless  velocity,  uj/u^l  +  x),  the  parameter 

X,  denotes  a  dimensionless  distance  related  to  the  x-  and  y-coordlnates  of  the 
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ellipse  by  x  »  1/2  ax  S,  y  ■  ax  S  measured  from  the  nose,  and  (=a/x) 
represents  a  dimensionless  angle  of  attack.  The  parameter  i  Is  also  related 
to  the  surface  distance  s  by 


Me  next  define  a  dimensionless  distance  n  by 


n 


r  R(1  +  t)  i^2  n 

SiT?i  * 


with  R  =  2a n»/v,  and  a  dimensionless  stream  function  f  by 

f(E,n,t)  =  [(1  +  T)aueovT2]’1^2^(s,n,t) 


(2.6a) 


(2.6b) 


Introducing  these  relations  together  with  Eq.  (2.5)  into  Eqs.  (2.1)  and  (2.2), 
it  can  be  shown  that  the  continuity  and  momentum  equations  can  be  written  as 


1  + 


.  1 2 


.2.— 


3tl  A  /  A  3U  .  A| 

+  "  t5‘>“e3rt('  +  5  )  lr 


f" 


3f 

3E 


+  (1  +  E2)|^-  (2.7) 
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Here  primes  denote  differentiation  with  respect  to  n  and  tj  *  (1  +  xJuJVaT  . 

The  boundary  conditions  for  f  and  f'  become 

f  ■  f'  *  0  at  n  =  0, 

f‘  -  (1  +  S2)1/2  UgtC.t,)  as  n-»  (2.8) 


The  definition  of  ue(s,tj)  Is  given  by 


ue(C,t) 


ue(E.t) 

Ujl  +  T) 


*  +  w 

/  1  +  S2 


(2.9) 


To  complete  the  formulation  of  the  problem,  initial  conditions  must  be  speci¬ 
fied  In  the  (t,n)  plane  at  some  s  =  sQ,  either  on  the  lower  or  upper  surface 
of  the  airfoil  as  well  as  Initial  conditions  In  the  plane  on  both  surfaces  of 
the  airfoil.  In  the  latter  case.  If  we  assume  that  steady-flow  conditions 
prevail  at  t  *  0,  then  the  Initial  conditions  In  the  (s,n)-plane  can  easily  be 
generated  for  both  surfaces  by  solving  the  governing  equations  for  steady  flow 
which  In  this  case,  are  given  by  Eq.  (2.1)  and  by 


There  is  no  problem  with  the  initial  conditions  for  Eqs.  (2.1)  and  (2.10)  since 
the  calculations  start  at  the  stagnation  point. 

The  generation  of  the  initial  conditions  in  the  (t.n)-plane  are  not  so 
straightforward  to  obtain  as  is  discussed  in  Section  3.1. 

3.0  SOLUTION  PROCEDURE 

The  solution  procedure  for  the  set  of  equations  and  boundary  and  initial  con¬ 
ditions  given  in  Section  2  can  be  achieved  in  two  parts  concerned,  respec¬ 
tively,  with  the  leading  edge  and  downstream  region.  These  two  parts  are 
considered  in  the  following  two  subsections.  In  both  cases  the  solution 
procedure  makes  use  of  Keller's  box  method,  which  is  a  two-point  finite- 
difference  scheme  extensively  used  for  the  solution  of  parabolic  partial- 
differential  equations,  as  discussed  by  Bradshaw  et  al.  . 

3.1  Leading-Edge  Region 

The  generation  of  the  initial  conditions  in  the  (t,n)-plane  at  s  =  sQ 
requires  a  special  numerical  procedure.  Given,  as  we  are,  the  complete  veloc¬ 
ity  profile  distribution  on  the  previous  time  line,  there  is,  In  principle,  no 
difficulty  In  computing  values  on  the  next  time  line  by  an  explicit  method,  but 
if  we  wish  to  avoid  the  stability  problems  associated  with  such  a  method  by 
using  an  implicit  method,  we  are  Immediately  faced  with  the  problem  of  gener¬ 
ating  a  starting  profile  on  the  new  time  line. 

In  order  to  explain  the  problem  further,  it  is  instructive  to  see  what  happens 

to  the  stagnation  point  as  a  function  of  time.  For  this  purpose  let  us 

consider  Eq.  (2.4).  If  we  choose  SQ(tj)  to  be  of  the  form  eq(1  +  A  slnutj), 

and  let  u°  =  u°/u  (1  +  t),  then  Eq.  (2.4)  becomes 
e  e  ® 


K  +  50U  +  A  slnut.) 
u^,tT) - 2 - - - L- 

✓  1  +  *z 

where  u  is  related  to  the  dimensional  frequency  by  ^  ~T  j  \  i 


(3.1) 


Here  E0  and  A  denote  parameters  that  need  to  be  specified.  Since  by 
definition  u£  *  0  at  the  stagnation  point.  Its  location,  Es,  is  given  by 

Cs--C0(1+A8lnai>  (3.2) 

and  so  the  upper  and  lower  surfaces  of  the  airfoil  as  functions  of  time  are 
defined,  in  particular,  by  £  >  and  K  <  E$.  For  example,  let  us  take  A  «  1, 
u  *  n/4  and  plot  5s/50  in  the  (t,E)-plane,  as  shown  in  Figure  1  for  one  cycle 
(0<  t,  <  8). 


fe 


1 


Fig.  1.  Variation  of  stagnation  point  with  time  for  one  cycle  according  to 
Eq.  (3.2),  with  u  ■  ir/4,  A  »  1. 

When  tj  =  2,  the  stagnation  point  Is  at  -2EQ,  whei  tj  «  6  Is  at  0,  etc.  If 
Ks  were  fixed,  we  could  assume  that  u  ■  0  at  E  *  s0  for  all  time  and  all  n> 
but  this  Is  not  the  case.  It  Is  also  possible  to  assume  that  the  stagnation 
point  Is  coincident  with  zero  u-veloclty  for  a  prescribed  time.  However,  we 
should  note  that  the  stagnation  point  given  by  Eq.  (3.2)  Is  based  on  the  van¬ 
ishing  of  the  external  velocity.  For  a  time- dependent  flow,  this  does  not 
necessarily  Imply  that  the  u-veloclty  Is  zero  across  the  layer  for  a  given 
E-locatlon  and  specified  time.  This  point  Is  substantiated  by  the  results 
shown  In  Figure  2  taken  from  Ref.  8  and  obtained  with  a  novel  numerical  pro¬ 
cedure  called  the  characteristic  box  scheme.  It  is  also  evident  from  Figure  2 
that  flow  reversals  do  occur  due  to  the  movement  of  the  locus  of  zero 
u-velocity  across  the  layer.  This  causes  numerical  instabilities  which  can  be 
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avoided  by  using  either  the  zig-zag  box  or  the  characteristic  box  finite- 
difference  scheme.  The  details  of  these  numerical  schemes  have  been  reported 
In  Ref.  7  and,  with  special  reference  to  oscillating  airfoils.  In  Ref.  8. 

3.2  Downstream  Region 

A  solution  to  the  leading-edge  region,  obtained  by  the  procedure  of  Section 
3.1,  may  be  used  as  Initial  conditions  for  the  solution  of  the  system  of  equa¬ 
tions  given  by  Eqs.  (2.7)  and  (2.8).  This  can  conveniently  be  done  by  using 

Keller's  box  scheme,  which  reduces  Eq.  (2.7)  to  a  first-order  system.  With  f' 
and  f"  represented  by  r  and  v,  respectively,  we  write 

f  =  r  (3.3a] 

r'  =  v  (3.3b! 


and  obtain 

v' +  TTp- ,!'1"  5z)«e  ♦  (i  ♦  c2>3/2 

-  >■  If  - v  (1  +  cZ)  %  (3-3c 

With  this  notation,  the  boundary  conditions  given  by  Eq.  (2.8)  can  be  written 
as 

f  ■  r  ■  0  at  n  *  0  (3.4) 

re  ■  (1  +  £Z)1/2ue  at  n  *  ne  (3.5) 

The  system  of  Eqs.  (3.3)  and  (3.4)  has  been  solved  by  the  numerical  procedure 
of  Ref.  7. 


4.0  NATURE  OF  THE  SINGULARITY  FOR  AN  OSCILLATING  AIRFOIL 

One  phase  of  the  calculations  for  the  oscillating  airfoil  was  carried  out  by 
choosing  »  1,  A  ■  -1/2  and  u  •  0.1.  With  these  choices  the  maximum  value 
of  aeff,  defined  by 
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(4.1) 


<*eff  *  5<>n  +  A  sinutj) 

Is  sufficient  to  provoke  separation  with  a  strong  singularity  If  the  boundary 
layer  were  steady.  At  present  50,  A,  u  are  being  varied  to  examine  their 
effect  on  the  nature  of  singularity. 

The  unsteady  flow  calculations  displayed  in  Fig.  3  show  that  the  boundary  layer 
eventually  separates,  the  flow  remaining  smooth.  However,  just  downstream  of 
separation,  it  is  evident  that  a  singularity  develops  in  the  solution  in  the 
neighborhood  of  C  *  2.12  and  uti  =  308.75°  and  that  it  is  not  possible  to 
continue  the  solution  beyond  this  time  without  conceptual  changes  in  the 
mathematical  and  physical  formulation  of  the  problem.  While  this  is  a  satis¬ 
fying  conclusion,  and  may  be  interpreted  as  giving  theoretical  support  to 
experimental  observations  of  dynamic  stall,  it  should  be  treated  with  some 
caution.  Boundary- layer  singularities  have  been  the  subject  of  much  contro¬ 
versy  in  recent  years  and  it  is  clearly  important  to  make  sure  that  any  irreg¬ 
ularities  in  a  computed  solution  are  not  creatures  of  the  numerical  method 
used.  We,  however,  feel  confident  that  the  calculations  reported  here  are 
accurate  and  that  the  singularity  is  real. 

Figure  3a  shows  that  the  variation  of  the  displacement  thickness 

**  *  ?  FT  <4-2> 

is  generally  smooth  except  in  the  neighborhood  of  %  »  2.12  and  for  u>t]  = 

308.75°.  The  first  sign  of  irregularity  is  the  steepening  of  the  slope  of 
when  ut]  ■  300°.  A  local  maximum  of  occurs  at  5  *  2.12  when  uti  = 

308.75°.  When  the  same  results  are  plotted  for  a  displacement  velocity,  (d/d£) 
(ue5*),  (Fig.  3b),  we  observe  that  the  steepening  of  the  displacement  velocity 
near  %  =  2.12  is  dramatic.  For  example,  for  uti  ■  300°,  the  peak  is  at  5  = 

2.125,  for  uti  *  305°,  it  is  at  5  ■  2.105,  for  wtj  «  307.5°,  it  is  at  5  « 

2.09  and  finally  for  ut]  *  308.75°,  the  peak  moves  to  X,  ■  2.08.  It  should  be 
noted  that  the  maximum  value  of  displacement  velocity  moves  towards  the  separation 
point  with  increasing  bit];  the  same  behavior  will  be  shown  to  occur  for 
the  circular  cylinder  discussed  below. 


f.-; 


As  shown  in  Fig.  3c,  the  wall  shear  parameter  f#  shows  no  signs  of  irregularity 

for  uti  <  308.75°  but  a  deep  minimum  In  f#  occurs  near  5  =  2.15,  i.e.  near 
the  peak  of  ft*. 

It  is  interesting  and  useful  to  compare  the  results  presented  in  Fig.  3  for  an 
oscillating  airfoil  with  those  obtained  for  a  circular  cylinder  started 
impulsively  from  rest.  This  comparison  lends  support  to  the  accuracy  of  the 
present  calculation  method  and  at  the  same  time  enables  us  to  compare  the 
characteristics  of  two  distinctly  different  unsteady  flows  near  the  singularity 
location.  The  circular  cylinder  problem  has  been  extensively  studied  as 
reported  in  Refs.  9  and  10  and  the  present  results  shown  below  are  in  close 
agreement  with  those  of  previous  authors,  but  with  subtle  differences  which  may 
have  important  implications. 

Figure  4  shows  the  results  obtained  by  Cebeci2  for  the  circular  cylinder 
problem.  As  in  the  case  of  the  oscillating  airfoil,  the  flow  separates  and 
remains  smooth  up  to  the  separation  point.  However,  just  downstream  of  sepa¬ 
ration  with  increasing  time,  a  singularity  develops  in  the  neighborhood  of  9 
=  112°  and  t  =  3.0  and  it  was  not  possible  to  continue  the  boundary- layer 
calculations  beyond  this  time  and  angular  location.  From  Fig.  4a  we  see  that 
while  the  variation  of  displacement  thickness  is  smooth  for  values  of  9  less 
than  108°,  it  begins  to  steepen  dramatically  thereafter.  The  same  results  are 
plotted  in  Fig.  4b  to  demonstrate  that,  as  in  Fig.  3b,  the  displacement  veloc¬ 
ity  exhibits  a  maximum  which  Increases  rapidly  with  time.  Again  the  maximum 
shifts  towards  the  location  of  separation  with  increasing  time. 

The  results  of  local  skin-friction  coefficient  calculations  in  Fig.  4c  follow 
similar  trends  of  those  obtained  for  the  oscillating  airfoil.  In  both  cases, 
the  distributions  pass  through  zero  with  no  signs  of  Irregularity  and  do  not 
exhibit  any  breakdown  before  the  time  corresponding  to  the  singularity. 


* 
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£  The  corresponding  displacement  thickness  results  of  Cebec1c  together  with  the 

new  calculations  of  the  velocity  profiles  are  reproduced  for  comparison 
H  purposes.  As  can  be  seen  from  Fig.  5,  the  agreement  between  the  sets  of  calc- 

J?  ulatlons  for  three  values  of  t  •  2,  2.5  and  2.75  Is  excellent.  The  velocity 


The  very  careful  calculations  of  van  Dommelen  and  Shen  are  reproduced  on 
Figures  5  and  6  for  displacement  thickness  and  velocity  profiles,  respectively. 
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profiles  of  Fig.  6a,  which  correspond  to  a  time  t  *  2.75  as  in  Fig.  5c,  are 
also  in  excellent  agreement  for  various  angular  locations.  In  contrast,  the 
calculated  velocity  profiles  of  van  Dommelen  and  Shen11  at  t  *  2.984375  show 
differences  from  the  present  results  obtained  at  t  =  2.9875.  The  figure  con¬ 
firms  the  expected  close  agreement  of  the  two  sets  of  results  at  the  two 
smallest  angular  locations,  but  indicates  significant  differences  at  the  two 
highest  values.  The  trend  is  different  in  that  the  present  results  show  that 
the  location  and  the  magnitude  of  the  maximum  negative  velocity  increases  with 
angular  location.  Also  the  tendency  for  flattening  of  the  velocity  profiles 
in  the  vicinity  of  the  singularity  is  not  confirmed  by  the  present  results. 

Figure  7  shows  the  velocity  profiles  obtained  by  Cebeci2  for  two  values  of  t 
as  a  function  of  angular  location.  It  is  clear  that  the  magnitude  of  negative 
velocity  increases  with  angular  location  and  suggests  that  as  the  singularity 
is  approached,  the  magnitude  of  the  negative  velocity  will  tend  to  infinity. 

3 

Figure  8  allows  comparison  of  the  displacement  velocities  obtained  by  Cowley 
and  by  the  present  method  for  four  values  of  time.  We  would  expect,  from  the 
previous  comparisons  that  the  two  sets  of  results  would  be  in  close  accord  at 
least  for  times  up  to  2.75.  The  figure  shows  the  expected  close  agreement 
until  the  maximum  value  is  approached.  The  discrepancies  apparent  at  higher 
values  of  e  cannot  readily  be  explained,  and  it  should  be  noted  that  the 
location  and  time  of  singularity  occur  at  different  values  of  0;  the  results 
of  Cowley  and  van  Dommelen  and  Shen  appear  to  agree  in  this  respect.  The 
reasons  for  these  discrepancies  are  presently  under  Investigation. 
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Velocity  profiles  in  the  immediate  neighborhood  of  the  stagnation  line  at  different  times  for 
u  =  it/4  and  A  =  1.  The  dashed  lines  indicate  the  locus  of  zero  u-velocity  across  the  layer. 


parameter 


Computed  results  of  Cebeci  for  the  circular  cylinder,  (a)  Displacement 
(b)  Displacement  velocity,  d/de(u  <S*).  (c)  Wall  shear  parameter,  f". 


Fig.  5.  Comparison  between  the  displacement  thickness  values  obtained  by 

van  Dommelen  and  Shen*  (circles)  and  by  Cebecl2  (solid  line)  for  the 
circular  cylinder,  (a)  t  *  2.0,  (b)  t  »  2.5,  (c)  t  »  2.75.  (x  Is  In 
radians.) 


1  van  Pommel en  and  Shen  (solid  lines) 
van  Pommelen  and  Shen)  and  t  =  2.9875 


Comparison  between  the  displacement  velocity  v< 
present  method  (dashed  lines)  for  the  circular 


